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Abstract 

The  Quraishi-Fahidy  method  [Can.  J.  Chem.  Eng.  59  (1981)  563]  was  employed  to  derive  characteristic  dimensionless  numbers  for 
the  membrane-electrolyte,  cathode  catalyst  layer  and  gas  diffuser,  respectively,  based  on  the  model  presented  by  Bernardi  and  Verbrugge 
for  polymer  electrolyte  fuel  cells  [AIChE  J.  37  (1991)  1151].  Monomial  correlations  among  dimensionless  numbers  were  developed  and 
tested  against  experimental  and  mathematical  modeling  results.  Dimensionless  numbers  comparing  the  bulk  and  surface-convective  ionic 
conductivities,  the  electric  and  viscous  forces  and  the  current  density  and  the  fixed  surface  charges,  were  employed  to  describe  the  membrane 
ohmic  drop  and  its  non-linear  dependence  on  current  density  due  to  membrane  dehydration.  The  analysis  of  the  catalyst  layer  yielded  electrode 
kinetic  equivalents  of  the  second  Damkohler  number  and  Thiele  modulus,  influencing  the  penetration  depth  of  the  oxygen  reduction  front  based 
on  the  pseudohomogeneous  film  model.  The  correlating  equations  for  the  catalyst  layer  could  describe  in  a  general  analytical  form,  all  the 
possible  electrode  polarization  scenarios  such  as  electrode  kinetic  control  coupled  or  not  with  ionic  and/or  oxygen  mass  transport  limitation. 
For  the  gas  diffusion-backing  layer  correlations  are  presented  in  terms  of  the  Nusselt  number  for  mass  transfer  in  electrochemical  systems. 
The  dimensionless  number-based  correlating  equations  for  the  membrane  electrode  assembly  (MEA)  could  provide  a  practical  approach  to 
quantify  single-cell  polarization  results  obtained  under  a  variety  of  experimental  conditions  and  to  implement  them  in  models  of  the  fuel  cell 
stack. 

©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

It  has  been  forecasted  that  polymer  electrolyte  fuel  cells 
(PEFC)  could  reach  widespread  commercialization  over  the 
next  one  or  two  decades  as  low  emission  power  sources 
for  a  variety  of  applications  ranging  from  electronic  de¬ 
vices  to  transportation.  The  large-scale  commercialization  of 
fuel  cells  hinges  on  both  socio-economic  ‘pull’  and  tech¬ 
nology  ‘push’.  Both  factors  are  strongly  dependent  on  the 
performance-cost  relationship.  Mathematical  modeling  has 
been  an  integral  part  of  fuel  cell  research  and  development, 
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initially  at  the  membrane-gas  diffusion  electrode  assembly 
(MEA)  level  (coupled  or  not  with  transport  phenomena  in 
the  flow  channels)  [1-12],  followed  up  by  attempts  to  develop 
stack  models  with  varying  degree  of  completeness  [13-16]. 

It  is  beyond  the  objective  of  this  section  to  give  a  compre¬ 
hensive  review  of  these  models.  However,  it  must  be  noted 
that  the  modeling  effort  coupled  with  experimental  validation 
provided  a  fairly  good  phenomenological  understanding  of 
the  physico-chemical  and  transport  phenomena  occurring  at 
the  MEA  level  and  the  projected  impact  on  single-cell  per¬ 
formance.  The  mechanistic  understanding  of  ion  and  water 
transport  in  the  polymer  membrane-electrolyte  together  with 
first-principle-based  analysis  and  the  practical  application  of 
these  models  for  membrane  performance  predictions  in  oper- 


106 


E.L.  Gy  enge  /  Journal  of  Power  Sources  152  (2005)  105-121 


Nomenclature 

ac  specific  surface  area  of  the  catalyst  layer 
(m2  m-3) 

AC  acid  capacity  of  the  membrane  (eq  kg-1) 
b  Tafel  slope  for  O2  reduction  (V  dec-1) 

BW  membrane  basis  weight  (kg  m-2) 

Cb  bulk  concentration  (mol  m-3) 

Cf  fixed  charge-site  concentration  of  the  mem¬ 
brane  (mol  m-3) 

Z?a-b  binary  gas  diffusion  coefficient  (m2  s-1) 

/)]!  H+  diffusion  coefficient  in  the  bulk  membrane- 
electrolyte  (m2  s-1) 

Do,  dissolved  O2  diffusion  coefficient  in  the  bulk 
membrane-electrolyte  (m2  s-1) 

Do2jC  effective  dissolved  O2  diffusion  coefficient  in 
the  membrane -electrolyte  phase  of  the  catalyst 
layer  (m2  s-1) 

Dan  the  second  Damkohler  number  (electrochemi¬ 
cal  definition) 

Eoc  open-circuit  cell  voltage  (V) 

standard  potential  for  O2  reduction  (V  versus 
SHE) 

EW  membrane  equivalent  weight  (kg  eq-1) 

/  electrochemical  factor  ( =aF/RT)  (V-1) 

F  Faraday’s  constant  (=96,487  C  mol^  or 
Ceq-1) 

H  Henry’s  constant  (Pam3  mol-1) 

i  superficial  current  density  in  the  membrane- 

electrolyte  phase  (A  m-2) 

j'l  mass  transfer  limiting  superficial  current  den¬ 
sity  (A  m-2) 

is  superficial  current  density  in  the  solid  phase 

(Am-2) 

it  total  (operating)  superficial  current  density 

(Am-2) 

k?  membrane  hydraulic  permeability  (m2) 

k&  electrokinetic  permeability  (m2) 

k°  standard  heterogeneous  rate  constant  (ms-1) 

l  characteristic  length  (m) 

—  number  of  moles  of  electrons  per  mole  of 

species  j  (moleiectron  mol-1) 

Nd  number  of  dimensionless  groups 

As  1  number  of  SI  units 

Nue  mass  transfer  Nusselt  number  (electrochemi¬ 

cal  definition) 

Pd  inlet  air  pressure  in  the  diffuser  (N  m-2) 

A p  hydraulic  pressure  difference  (N  m-2) 

qf  fixed-negative  surface  charge  density  (C  m-2) 

rp  membrane  pore  radius  (m) 

R  universal  gas  constant  (J  mol-1  K-1) 

Re  Reynolds  number 

Sj  stoichiometric  coefficient  of  species  j  (either 

H+  or  02) 


Sh  Sherwood  number 

T  temperature  (K) 

u  ionic  mobility  (m2  s-1  V-1) 

i>0  superficial  water  velocity  (ms-1) 

Wa  Wagner  number 

x  gas  mole  fraction 

y  exponents  and/or  constants  in  the  correlating 

functions 

Zf  charge  number  for  the  fixed  sites  of  the  proton 

exchange  membrane  (=—  1) 

Z  position  along  the  MEA  thickness 

Zc  catalyst  layer  thickness  (m) 

Zd  diffuser  thickness  (m) 

Zm  ‘wet’  membrane  thickness  (m) 

Greek  symbols 

a  transfer  coefficient  for  O2  reduction 

yw  water  activity  coefficient 

i]  overpotential  for  O2  reduction  (V) 

9m  membrane  water  uptake  (m3  m-3) 

k  membrane-electrolyte  specific  conductivity 

(Sm-1) 

Ko  membrane-electrolyte  specific  conductivity  at 

zero  current  (i.e.  open-circuit  conditions) 

(Sm-1) 

X  membrane  water  content  (i.e.  number  of  water 

molecules  per  fixed  ionic  group) 

/x  dynamic  viscosity  (Pa  s) 

Vp  reaction  front  penetration  depth  inside  the  cat¬ 
alyst  layer  (m) 

7T  dimensionless  number 

p  density  (kg  m-3) 

a  electronic  conductivity  of  the  solid-matrix 

phase  (S  m-1) 

r  slope  of  the  <pE  versus  7rio  dependence 

(p  Thiele  modulus 

0  the  potential  in  the  ionic  conductor  polymer 

electrolyte  phase  (V) 

0s  the  potential  in  the  solid,  electronic  conductor, 
phase  (V) 

0  correlating  function 

Subscripts 

c  catalyst  layer  (it  also  refers  to  effective  values 

of  parameters  in  the  porous  catalyst  layer) 
d  gas  diffusion  layer  (it  also  refers  to  effec¬ 

tive  values  of  parameters  in  the  porous  gas- 
diffusion  layer) 

E  related  to  electrodes  and  electrochemical  sys¬ 

tems 

f  fixed-negative  charges  on  the  membrane 

1  limiting  value  for  linear  dependence  of  the 

membrane  ohmic  drop  versus  current  density 
L  mass  transfer  limiting 
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m  proton  exchange  membrane  (polymer  elec¬ 

trolyte) 

s  solid-matrix  electronic  conductor  phase 

t  total  (i.e.  electronic  and  ionic  conductor 

phases) 

w  water  phase 

0  inlet  to  the  gas-diffusion  layer 

1  catalyst/gas-diffusion  layer  interface 

2  bulk  membrane-electrolyte/catalyst  layer  in¬ 
terface 


ating  fuel  cells,  is  less  developed  and  it  is  a  subject  of  intense 
research  (for  a  review  see  [17]). 

With  the  fuel  cell  development  effort  focusing  increas¬ 
ingly  on  stack  engineering  and  optimization,  there  is  a  need, 
as  noted  also  by  Appleby  and  co-workers  [13],  for  com¬ 
prehensive  yet  numerically  less  extensive  equations  descrip¬ 
tive  of  the  MEA  performance  in  order  to  be  easily  imple¬ 
mented  into  complex  numerical  models  of  the  fuel  cell  stack. 
Therefore,  in  parallel  with  phenomenological  and  mechanis¬ 
tic  models,  empirical  equations  have  been  developed  to  de¬ 
scribe  the  fuel  cell  voltage  £  as  a  function  of  current  density  i. 
The  equation  proposed  by  Chamberlain  and  co-workers  (Eq. 
(1))  showed  good  fit  with  experimentally  measured  data  [18]: 

E  —  Eoc  —  b\ogi  —  Ri  —  yexp(a>0,  (1) 

where  b  is  the  Tafel  slope  for  oxygen  reduction  (V  dec-1), 
E  and  Eoc  the  operating  and  open-circuit  cell  voltages,  re¬ 
spectively  (V),  i  the  superficial  current  density  (Am-2),  R 
the  area-specific  ohmic  resistance  of  the  membrane  (Q  m2), 
y  (V)  and  <o  (m2  A-1)  the  empirical  curve  fitting  parameters. 

The  exponential  term  in  Eq.  (1)  expresses  empirically  the 
mass  transfer  limitation  of  the  polarization  curve  [18].  More¬ 
over,  the  coefficient  b  in  Eq.  (1)  is  often  used  as  an  additional 
curve  fitting  parameter. 

Amphlett  et  al.  developed  an  empirical  model  as  well,  to  fit 
the  experimentally  obtained  polarization  curves  of  a  Ballard 
fuel  cell  [19].  The  features  of  their  cell  voltage  model  were: 
linear  dependence  of  the  internal  resistance  on  current  density 
and  temperature,  and  logarithmic  dependence  of  the  overpo¬ 
tential  on  both  current  density  (i.e.  Tafel  regime)  and  effec¬ 
tive  oxygen  concentration  coupled  with  a  linear  temperature 
function  correlated  by  respective  parametric  coefficients. 

While  the  proposed  empirical  equations  are  successful 
in  describing  the  overall  single-cell  polarization  there  is  an 
inherent  lack  of  transparency  associated  with  them,  which 
could  restrict  the  applicability  to  the  specific  conditions  and 
possibly  cell  designs  used  in  determining  the  empirical  pa¬ 
rameters. 

The  goal  of  the  present  investigation  was  to  determine 
characteristic  dimensionless  numbers  for  each  component  of 
the  MEA  (exemplified  on  the  cathode  side.  Fig.  1)  and  fur¬ 
thermore,  to  evaluate  whether  monomial  correlations  among 


dimensionless  ratios  could  be  derived  to  provide,  after  testing 
against  experimental  data,  easy  to  use  equations  linking  es¬ 
sential  physico-chemical  properties  with  fuel  cell  operating 
parameters. 

Dimensional  analysis,  although  a  time-honored  technique 
(for  a  recent  review  see  the  monograph  by  Szirtes  [20]), 
as  a  stand-alone  modeling  method  for  chemical  engineer¬ 
ing  is  typically  limited  to  simple  cases  of  transport  processes 
where  a  list  of  variables  could  be  straightforwardly  compiled. 
This  approach  is  clearly  not  applicable  to  more  complex  sys¬ 
tems  especially  those  involving  both  transport  phenomena 
and  chemical/electrochemical  reactions.  However,  as  noted 
by  Churchill,  dimensional  analysis  applied  to  a  mathemati¬ 
cal  model  of  the  system,  eliminates  the  uncertainty  associated 
with  the  selection  of  variables  and  useful  quantitative  results 
in  the  form  of  dimensionless  correlations  could  be  obtained. 
In  those  cases  the  dimensional  analysis  results  could  describe 
the  system  to  the  extent  of  the  mathematical  model  they  are 
based  upon  [21]. 

In  a  review  of  dimensional  analysis  applications  for  elec¬ 
trochemical  reactors,  Selman  also  concluded  that  due  to  the 
large  number  of  interacting  variables,  ‘blind’  dimensional 
analysis  must  be  avoided  and  dimensional  analysis  could  only 
be  useful  when  coupled  with  a  mathematical  model  of  the  sys¬ 
tem  [22].  In  the  present  work  the  dimensional  analysis  pro¬ 
cedure  developed  by  Quraishi  and  Fahidy  [23]  was  applied 
to  the  variables,  parameters  and  constants  encountered  in 
the  one-dimensional  MEA  model  presented  by  Bernardi  and 
Verbrugge  [4,5].  Furthermore,  monomial  correlating  func¬ 
tions  among  dimensionless  numbers  are  proposed  and  tested 
against  experimental  data  and  mathematical  modeling  results 
from  the  literature. 


2.  Model  summary  of  the  membrane-gas  diffusion 
electrode  assembly 

Fig.  1  shows  schematically  the  three  components  of  the 
MEA,  gas  diffusion  layer,  catalyst  layer  (i.e.  reaction  zone) 
and  the  membrane-electrolyte.  Here  the  MEA  is  described  in 
relation  to  the  cathode  reaction: 

02  +  4H+  +  4e-  2H20, 

E°  =  1.23  VsHEat  298  K.  (2) 

The  Bernardi- Verbrugge  model,  selected  as  basis  for  the 
dimensionless  analysis,  has  been  discussed  extensively  in 
the  literature  and  compared  to  experimental  results  [4,5]. 
It  is  essentially  based  on  the  macroscopic,  volume  aver¬ 
aged,  porous  electrode  model  [24],  coupled  on  the  mem¬ 
brane  side  with  the  dilute  electrolyte  assumption  for  ion  trans¬ 
port  and  Schlogl’s  equation  for  water  transport,  whilst  the 
Stefan-Maxwell  equation  and  material  balance  equations  are 
employed  to  describe  the  gas  diffuser.  Ohm’s  law,  written  in 
terms  of  effective  conductivities,  is  applied  to  describe  the 
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from  flow  channel 


Unconsumed  gases;  H20 
to  flow  channel 
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Polymer  Catalyst  Gas 

electrolyte  layer  diffuser 


Fig.  1 .  Schematic  representation  of  the  membrane-gas  diffusion  electrode  assembly  (MEA)  on  the  cathode  side  of  the  polymer  electrolyte  fuel  cell. 


potential  drop  in  the  electronic  and  ionic  conductive  phases 
of  the  porous  electrode  as  a  function  of  current  density. 

The  transfer  current  per  unit  electrode  volume  between 
the  solid-matrix  and  polymer  electrolyte  phases  for  the  ca¬ 
thodic  direction  of  oxygen  reduction  (Eq.  (2)),  considering 
first-order  kinetics  with  respect  to  dissolved  oxygen  concen¬ 
tration  and  the  proton  concentration  effect  factored  in  the  rate 
constant,  is  given  by 

%  =  —Fack°Co2  exp (-/(<£s  -  0  -  E0)) 
d  Z  so2 

=  aci0  ( exp (-/(<2>s  -0-  E0)),  (3) 

\ko2,ref/ 

where  f  =  pfp,  ac  the  specific  surface  area  of  the  cata¬ 
lyst  layer,  Co2  the  actual  dissolved  O2  concentration  in  the 
polymer  electrolyte  phase  of  the  catalyst  layer,  Co2,ref  a 
reference  O2  concentration  at  which  the  exchange  current 
density  is  specified,  £°  the  standard  potential  for  the  O2 
electrode,  F  the  Faraday’s  constant,  k°  the  standard  het¬ 
erogeneous  rate  constant,  ;'o  the  exchange  current  density 
(70  =  (n /sO2)Fk0CO2iief),  n  the  number  of  electrons  involved 
in  O2  reduction,  so2  the  stoichiometric  coefficient  for  O2,  R 
the  universal  gas  constant,  T  the  temperature,  Z  the  position  in 
the  catalyst  layer,  a  the  cathodic  transfer  coefficient  for  O2  re¬ 
duction,  while  0S  and  0  the  potentials  in  the  solid  electronic 
conductor  and  polymer  electrolyte  phases,  respectively. 

Schlogl’s  equation  for  water  transport  in  the  membrane 
considers  water  movement  due  to  both  hydraulic  pressure 
gradient  dp/dZ,  and  electroosmotic  drag,  which  is  propor¬ 
tional  to  the  potential  gradient  across  the  membrane  d0/dZ 
[4]; 


Vm 


kp  dp  kcp  d0 

H - 

p  dZ  p  dZ 


(4) 


with  Cf  fixed  surface  charge-site  (i.e.  sulfonate  group)  con¬ 
centration;  kp,  kcp  the  hydraulic  and  electrokinetic  membrane 
permeability,  respectively;  Zf  the  charge  number;  um  the  net 


water  velocity  in  the  membrane  and  p  the  water  (i.e.  pore- 
fluid)  dynamic  viscosity. 

The  complete  set  of  model  equations  for  membrane,  cat¬ 
alyst  layer  and  gas  diffuser  are  presented  in  the  work  of 
Bernardi  and  Verbrugge  [4,5].  The  associated  list  of  vari¬ 
ables,  parameters  and  constants  as  applied  in  the  dimensional 
analysis,  are  presented  in  Tables  1-3. 


3.  Dimensional  analysis  technique 

In  the  present  work  the  Quraishi-Fahidy  method  [23]  has 
been  applied  to  identify  the  dimensionless  numbers  describ¬ 
ing  the  behaviour  of  the  proton  exchange  membrane,  catalyst 
layer  and  gas-diffusion  layer  according  to  the  model  outlined 
in  the  previous  section.  The  Quraishi-Fahidy  dimensional 
analysis  technique  comprises  the  following  steps:  (1)  a  list 
containing  all  the  variables,  parameters  and  constants  with 
their  respective  SI  units  (using  unity  for  dimensionless  quan¬ 
tities)  is  complied,  based  on  the  governing  model  equations 
(algebraic,  differential  and/or  integral);  (2)  each  dimensional 
quantity  on  the  list  is  divided  by  its  SI  units,  generating  quan¬ 
tity  per  unit  ratios  (or  groups);  (3)  considering  /Vs [  the  total 
number  of  independent  SI  units  involved  in  the  model,  a  num¬ 
ber  of  /Vs  1  quantity-unit  groups  are  selected  such  that  to  have 
the  simplest  possible  mathematical  form  and  all  SI  units  en¬ 
countered  to  be  represented  in  the  /Vsi  groups;  (4)  each  of 
the  Nsi  number  of  groups  is  assigned  a  value  of  1,  desig¬ 
nated  as  ‘primary’  groups;  (5)  the  rest  of  the  quantity-unit 
groups  on  the  list  are  converted  into  an  initial  set  of  dimen¬ 
sionless  numbers  jr by  substituting  the  units  using  the  Nsi 
primary  groups;  (6)  the  initial  set  of  dimensionless  numbers 
jr,  is  further  transformed  by  algebraic  manipulation  in  order  to 
obtain  as  many  as  possible  conventionally  known  dimension¬ 
less  numbers  and/or  to  generate  new  dimensionless  numbers 
with  clear  physico-chemical  meaning. 

Tables  1-3  illustrate  the  outlined  procedure  and  summa¬ 
rize  the  set  of  dimensionless  numbers  that  are  characteris- 
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Table  1 


Dimensional  analysis  of  the  membrane 


Number 

Quantity 

Primary  SI  units 

Quantity/units 

Dimensionless  numbers  ni 

Initial  jt/  Transform 

Final  7T; 

i 

i  Nt) 

Am-2 

.  2 

i  m 

A 

i 

- 

2 

Zm 

m 

Zm 

m 

i 

- 

3 

\zf\Ct 

mol  m~3 

IZflCfin3 

mol 

i 

- 

4 

P 

kgm-3 

pm3 

kg 

i 

- 

5 

P 

kgm-1  s_1 

pms 

i 

- 

6 

Da+ 

2  -1 
m  s 

dh+  s 
m2 

dh+p  X 

P  JT1 

711  = 

7 

(*)f 

As  mol-1 

nFmol 

■SH+ 

»F|zf|CfM,  1 

iZmsn+  P  ^2 

iZms  h+  P 
-  nF\zf\Cftl 

8 

Co2 

mol  mT3 

Co,  m3 

mol 

Co2 

liflCf 

Co2 
~  |Zf|Cf 

9 

kpAp 

kg  m  s-2 

(kpAp)s2 

kgm 

(kpAp)p 

P2 

_  (keAp)p 

7T4  —  2 

p2 

10 

k^Acp 

kgm4  A-1  s-3 

(k<pA<P)  As3 
kgm4 

(k@  A<P)ip2  Zm 

^3  7T57T27T4 

Tr  _  n\zf\CfF(k<pA<P) 

5  sH+kpAp 

11 

12 

K 

Do, 

A2  s3  kg-1  m~3 

m2  s_1 

/ckgm3 

KP?  7T6 

A2  s3 

Do2  s 

i2Z4mP 2  7T2 

do2p  1 

p  7T7 

T  -  _ 

6  (,,rl|zf|cf)2z^1 

tic  for  the  membrane. 

catalyst  and  gas-diffusion  layer,  re- 

alyst  layer  by  1 1  numbers,  jtq,  : r§, . 

. .,  7Ti7,  while  the  gas- 

spectively.  For  all  three  cases  Ns \  =  5  (i.e.,  kg,  m. 

,  s,  mol 

diffusion  (or  backing)  layer  could  be  described  using  8  num¬ 

and  A). Thus, 

the  proton  exchange  membrane  could  be  char- 

bers,  7ri8, . . .,  7T25.  The  physico-chemical  significance  and 

acterized  by  7  dimensionless  numbers,  n\, . .  .,777, 

the  cat- 

implications  of  the  obtained  dimensionless  numbers  for  the 

Table  2 

Dimensional  analysis  of  the  cathode  catalyst  layer 

Number 

Quantity 

Primary  SI  units 

Quantity/units 

Dimensionless  numbers  TTj 

Initial  ici 

Transform 

Final 

i 

i 

Am'2 

.  2 

i  m 

A 

i 

- 

- 

2 

Zc 

m 

Zc 

m 

i 

- 

- 

3 

Co 2 

mol  m-3 

Co2  m3 
mol 

i 

- 

- 

4 

P 

kgm-3 

pm3 

kg 

i 

- 

- 

5 

P 

kg  m  1  s_  1 

pms 

kg 

i 

- 

- 

6 

e-/(®s-#-E°) 

- 

e— /(#S -*-E°) 

- 

7ZQ  = 

e-/,fs -4>-E°) 

7 

Dq2,c 

2  —1 
m  s 

Do2,c  s 

D02,c  P 

Jig  = 

P 

m2 

P 

*8 

Do2,cP 

8 

ack° 

s_  1 

( ack° ) s 

(ack°)pZl 

JT9 

JCg  = 

(ac  k°)Zl 

P 

*8 

D02,c 

9 

(jL. \F 

As  mol-1 

nFmol 

hFCq2p 

_ZE8_ 

7T10  = 

iZcso2 

V°2) 

jo2  As 

so2iZcp 

^10 

hFDq2  cCq2 

10 

h 

Am~2 

.  2 

1 1  m 

A 

h 

i 

- 

7Tn  = 

:  It 
i 

11 

\zt\Cf 

mol  m-3 

|Zf|Cfm3 

mol 

kflCr 

Co2 

- 

H 12  = 

.  IztlCf 

Co2 

12 

kp.c  Ap 

kg  m  s-2 

( kp,cApc)s 2 
kgm 

( ke,cApc)p 

M2 

- 

H 13  = 

(kp,cApc)p 

P2 

13 

fcpjCA# 

kgm4  A-1s-3 

(r®,0  A#)  As3 
kgm4 

(r®cA<f>)ip2zc 

f3 

77147713^10 

^8*12 

7Ti4  = 

n  |zf  |  CfF(k$tC  A4>) 
sO,  kps  Ap 

14 

A0S 

kgm2  A-1  s-3 

A$sAs3 

A<PsiZ2p2 

771577137710 

7T15  = 

n|zt|CfFA(fsZ2 

kg  m2 

P? 

^8^12 

so2kp,cA.p 

15 

Kc 

A2  s3  kg-1  m-3 

kc  kg  m3 

KC  P 3 

7Ti6  = 

Kcps2^ 

A2  s3 

i2Z4p2 

'  (n|zf|FCf)2Z2 

16 

A2  s3  kg-1  m-3 

ac  kg  m3 

p3 

M’S)' 

7Z\1  = 

Pcpsl2 

A2  s3 

i2Z^p2 

'  (n|zf|FCf)2Z2 
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Table  3 

Dimensional  analysis  of  the  gas  diffusion-backing  layer 

Number  Quantity  Primary  SI  units  Quantity/units  Dimensionless  numbers  m 


1 

h  (=*t) 

Am-2 

2 

is  m 

'a 

2 

Zd 

m 

Zd 

m 

3 

P 

kg  m-3 

pm3 

kg 

4 

p 

kg  m-1  s-1 

pms 

kg 

5 

0 

X 

< 

mol  m-3 

Axq2  Pd  m3 
RT  mol 

6 

Dti2-o2,i 

2  —l 
m  s 

£)N2— 02,d  s 
m2 

7 

Dw— air,d 

2  —l 
m  s 

77 w — air.d  S 
m2 

8 

v  (_  Pw,sat\ 

*w-sat  1“  Pd  ) 

- 

- 

9 

*n2 

- 

- 

10 

vo,d 

ms-1 

^0,ds 

m 

11 

(*)f 

As  mol-1 

nFmol 

5q2  As 

12 

A&s 

kgm2  A^1  s-3 

A0S  As3 
kgm2 

13 

A2  s3  kg-1  m~3 

<Tdkgm3 

A2s3 

Initial  m  Transform  Final  7T( 


1 

1 

1 

1 

1 


°n2-o2  'AP  l 

F  *18 

^vv— air.dP  1 

P-  7TJ9 


•*w,sat 


^18 

7Ti9 


P 

^2-02,  AP 
V- 

£>w— air,d  P 


1*20  —  -*\v,sat 


*N2 

v0,dPZd 

P 

nF Axq2  pdp 
isZdso2  pRT 
Atf>sisZ3p2 

°dP3 

t2zy 


^18 

^23 


.  ^19 
^23,w 


^24^22^23 

^18 


*25 -T 


K  21  : 
7T22  : 


■  xn2 

vQ,dPZd 

P 


hso7  ZdRT 

7T93  =  rr. - t -  Or  71t\  \ 

nFD^2-02,dAx02Pd  ’ 


isso2  ZdRT 

nFDw_air2,dAxo2  pd 


^24  = 
^25  = 


”FA*^  (Axo, 
^o2uo,d^  V  W2 


Pd_\ 
RT  ) 


sn.  /7Aa'°2  flr 


membrane-gas  diffusion  electrode  assembly  analysis  are  dis¬ 
cussed  in  the  next  paragraph. 

Furthermore,  the  relation  among  dimensionless  numbers 
(i.e.  correlating  equations)  can  be  represented  as 

n„  =  V{7tj \j  =  1, (Nd  -  1)},  n^j,  (5) 

where  /Vj  is  the  number  of  dimensionless  groups  describing 
the  system,  jrn  a  dimensionless  number  containing  the  de¬ 
pendent  variable,  and  'F  a  correlating  function,  that  is  either 
monomial  or  non-monomial  [20].  The  monomial  (i.e.  power 
series)  form  is 

=  yon{l7t22  ■  ■  ■  kna- iWd_1  (6) 

with  y  j  \j  —  1, . . . ,  (N&  —  1)  numeric  exponents  and  yo  a 
constant  coefficient.  The  numeric  exponents  and  constant 
could  be  determined  by  regression  using  experimental  data 
relating  the  dependent  variable  entering  the  expression  of  nn 
to  independent  variables  and  parameters  composing  tcj,  cou¬ 
pled  with  asymptotic  and  speculative  analysis  as  shown  by 
Churchill  and  co-workers  [25-27]. 

The  case  when  is  of  non-monomial  form  (i.e.,  it  in¬ 
cludes  one  or  more  of  the  following:  addition,  subtraction,  or 
transcendental  functions  of  nf)  leads  to  multiple  coefficients 
yo  (at  least  two  [20]).  Therefore,  the  quantitative  interpreta¬ 
tion  and  analysis  of  the  dimensionless  correlation  becomes 
more  difficult  as  compared  to  the  monomial  function  situa¬ 
tion.  To  detect  a  non-monomial  'F,  the  typical  procedure  is 
to  assume  initially  a  monomial  form  and  to  prove  it  leads  to 
a  contradiction  or  impossibility  [20]. 


In  the  next  paragraph  the  characteristic  dimensionless 
numbers  and  the  applicability  of  monomial  correlating  equa¬ 
tions  for  the  membrane,  catalyst  layer  and  gas-diffusion  layer, 
are  discussed. 


4.  Results  and  discussion 

4.1.  Dimensionless  numbers  and  monomial  correlating 
equations  for  the  membrane 

Table  1  shows  the  seven  dimensionless  numbers,  which 
are  characteristic  for  the  membrane,  obtained  by  applying 
the  Quraishi-Fahidy  technique  to  the  variables  and  parame¬ 
ters  encountered  in  the  model.  The  dimensionless  ratios  jti 
and  7t~]  are  membrane  Schmidt  numbers  relating  the  pore- 
fluid  (i.e.  water)  viscosity  to  the  proton  and  dissolved  O2 
diffusion  coefficients,  respectively,  whilst  713  is  simply  the 
concentration  ratio  of  dissolved  O2  and  fixed  charges. 

The  dimensionless  number  ni  (Table  1  and  Eq.  (7))  on 
the  other  hand  is  new,  and  it  relates  the  current  density  i, 
membrane  thickness  Zm->  and  the  fixed  surface  charge  con¬ 
centration  \zf\Cf  (which,  based  on  electroneutrality  is  equal 
to  the  membrane  proton  concentration,  |zf|Cf  =  CH+),  with 
the  pore-fluid  dynamic  viscosity  p  and  density  p.  Moreover, 
TT2  could  be  expressed  as 

<Z,n.Y|  [- p  i  Zn)S]  |  f  Zfjj+p  Jtj 

Jt2  =  -  = - =  — —  .  (/) 

nF\zt\Cfp  nFDli+Cii+  p  7Ti 
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In  Eq.  (7),  the  ratio  n*  has  an  expression  mathematically 
similar  but  conceptually  different  than  the  mass  transfer  Nus- 
selt  number  for  electrochemical  systems  Nuv,  introduced  by 
Ibl  [28]: 


The  dimensionless  ratio  jt g  on  the  other  hand,  relates  the 
bulk  (or  actual)  specific  ionic  conductivity  of  the  membrane 
k  to  the  surface-convective  conductivity,  which  is  related  to 
the  fixed  surface  charge  density  qf  (Cm-2): 


Nu  E 


il 

nFDAC’ 


(8) 


where  AC  is  the  concentration  difference  between  the  bulk 
and  electrode  interface,  D  the  diffusion  coefficient,  I  the  char¬ 
acteristic  length  and  n  the  number  of  electrons.  For  mass 
transfer  limiting  conditions,  i.e.  AC=Cb  and  /=;'l,  the  di¬ 
mensionless  ratio  given  by  Eq.  (8)  is  commonly  referred  to 
as  the  Sherwood  number  [23,29]. 

Considering  that  the  proton  concentration  in  the  mem¬ 
brane  is  constant,  the  interpretation  given  for  Nue  is  not 
applicable  to  tc%-  Substituting  the  Nernst-Einstein  equation, 
Z)H+  =  “h  f  RT ,  which  is  based  on  the  commonly  employed 
infinite-dilution  assumption,  shows  that  jt*  expresses  the  ra¬ 
tio  between  the  total  (i.e.  migration  plus  convection)  proton 
flux  in  the  membrane  pore  fluid  j  and  the  migra¬ 

tion  flux  corresponding  to  a  unitary,  1  V,  potential  difference 
across  the  membrane 

It  must  be  noted  that  DH+  increases  with  membrane  wa¬ 
ter  content  [30],  which  in  turn  is  a  function  of  the  operating 
conditions  of  the  fuel  cell.  Furthermore,  the  water  content 
across  the  membrane  thickness  might  not  be  uniform  [31] 
(e.g.,  the  cathode  side  could  be  more  hydrated  than  the  an¬ 
ode  side)  hence,  Z)H+  and  implicitly  7r|  could  have  different 
values  across  the  thickness  of  the  membrane.  The  dimension¬ 
less  number  jtj  on  the  other  hand  (Eq.  (7)),  is  constant  under 
given  operating  conditions. 

The  dimensionless  number  714  in  Table  1  is  a  ratio  of  iner¬ 
tial  to  viscous  forces  with  respect  to  the  membrane  pore- 
fluid,  where  the  inertial  force  is  related  to  the  hydraulic 
permeability-pressure  difference  product  kpAp  (Darcy’s 
law).  Thus,  714  is  a  membrane  pore-fluid  Reynolds  number, 
Rem- 

Additionally,  two  new  dimensionless  numbers  have  been 
identified,  which  were  not  discussed  previously  in  the  lit¬ 
erature,  JT5  and  7T6,  respectively  (Table  1).  The  expression 
of  7T5  includes  the  potential  difference  (or  drop)  across  the 
membrane-electrolyte  A  <t>,  in  relation  with  certain  physico¬ 
chemical  properties  of  the  membrane  such  as  electrokinetic 
and  hydraulic  permeabilities,  k@  and  kp,  and  fixed-surface 
charge  concentration  |zf|Cf.  Since  /;  =  4  and  .sH+  —4  (Eq. 
(2)),  7 r5  is  given  by 


_  "h+  ch+  ) 
Z, 


\zf\FCf(k0A‘t>)  electrical  force 

7T5  =  - } — 7 -  =  — - 7T - ■  (9) 

kpAp  inertial  force 

Thus,  715  could  be  interpreted  as  the  ratio  between  the 
pore  water  velocity  induced  by  the  electrical  force  (electro- 
osmotic  drag)  and  the  velocity  due  purely  to  the  inertial  force 
(see  also  Eq.  (4)).  A p  is  defined  here  as  the  hydraulic  pres¬ 
sure  difference  between  the  cathode  and  anode  side  of  the 
membrane. 


Kp  Kp 

(| Zf\FC{)2Zi  =  ^ 

bulk  (actual)  proton  conductivity 
surface-convective  proton  conductivity 


(10) 


In  the  analysis  of  electrokinetic  phenomena,  Newman  in¬ 
troduced  a  dimensionless  ratio  conceptually  similar  to 

a2 

i.e.,  where  q{\\  is  the  surface  charge  density  in  the  diffuse 
double  layer  and  k  is  the  specific  conductivity  at  the  center  of 

the  capillary.  The  surface-convective  conductivity  is  propor- 

2 

tional  to  f  and  it  is  due  to  the  convective  transport  of  ions 
by  the  fluid  moving  in  the  diffuse  layer  under  the  action  of 
the  electric  field  [32].  This  concept  is  applicable  to  the  pro¬ 
ton  exchange  membrane,  where  the  convective  surface  proton 
conductivity  can  be  considered  proportional  to  the  square  of 
the  fixed  surface  charge  density  qf. 

Thus,  7T6  compares  the  proton  conductivity  of  the  mem¬ 
brane  k,  which  is  dependent  on  the  actual  operating  con¬ 
ditions  (e.g.,  membrane  water  content  and  temperature  as 
shown  by  Eq.  (11)),  with  the  intrinsic  surface-convective 
conductivity  of  the  membrane  determined  by  the  physico¬ 
chemical  properties  of  the  membrane. 

Based  on  the  studies  of  Zawodinski  et  al.  [33],  Meier  and 
Eigenberger  [31]  established  the  following  equation  describ¬ 
ing  the  water  content  X  and  temperature  T  dependence  of  k 
(Sm-1): 


k  —  (0.46k  —  0.25)  exp 


1 


298.15 


0  <  X  <  30. 


(11) 


However,  in  operating  fuel  cells  it  is  difficult  to  determine 
experimentally  either  the  in  situ  water  content  or  the  specific 
conductivity  of  the  membrane.  Furthermore,  the  membrane 
specific  conductivity  could  be  anisotropic  across  its  thickness 
[31,34].  Therefore,  in  the  absence  of  relevant  data  under  fuel 
cell  operating  conditions,  in  the  expression  of  jTq,  k  could  be 
approximated  with  the  in  situ  specific  conductivity  at  zero 
current  (i.e.,  at  open-circuit  conditions)  kq. 

Biichi  and  Scherer  reported  for  Nafion®  membranes  of 
llOOgeq-1  equivalent  weight,  that  the  in  situ  specific  con¬ 
ductivity  at  zero  current  kq,  was  approximately  constant  at 
given  temperature  and  humidification  conditions  for  mem¬ 
brane  thickness  greater  than  about  150 |xm  (e.g.,  10 Sm-1 
at  333  K,  dry  O2  feed/humid  H2,  for  N105  and  N117  mem¬ 
branes  of  ‘wet’  thickness  163  and  203  pm,  respectively)  [34]. 
However,  the  same  authors  observed  that  for  thinner  mem¬ 
branes  (i.e.,  ‘wet’  thickness  less  than  about  150  pm  such  as 
N1 12:  60  pm,  NE  1135:  88  pm  and  N1 15:  148  pm)  the  thin- 
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ner  the  membrane  the  lower  the  specific  conductivity  at  zero 
current  [34] .  Discounting  potential  experimental  errors,  these 
differences  were  explained  by  conductivity  anisotropicity  at 
the  membrane  boundary  layer  adjacent  to  the  gas-diffusion 
electrode  structure  [34]. 

The  seven  dimensionless  numbers  given  in  Table  1  in  con¬ 
junction  with  a  correlating  function  could  be  employed  to 
describe  and  analyze  the  membrane  behaviour.  To  simplify 
the  correlating  function,  the  effect  of  the  dissolved  O2  con¬ 
centration  on  the  membrane  potential  drop,  expressed  by  7T3 
and  jtj,  was  neglected.  This  assumption  implies  a  complete 
consumption  of  O2  in  the  catalyst  layer. 

Employing  a  monomial  correlating  function  (Eq.  (6)),  and 
considering  A0  the  dependent  variable,  the  membrane  be¬ 
haviour  could  be  described  as 


^5  =  (12) 

(\zi\FCf(k0A&)\  _  (  iZmsH+p  Y2 

V  kvAp  )  -v°'m  |Zf|  Cffi) 

<*> 

The  correlating  equation  (13)  describes  in  a  dimensionless 
form,  the  potential  drop  across  the  membrane  as  a  function 
of  the  open-circuit  specific  conductivity,  membrane  physico¬ 
chemical  properties  and  fuel  cell  operating  conditions.  Next 
examples  are  given  for  the  determination  of  yi,  >’4,  yg  and 
yo.m  using  literature  studies  of  the  membrane  ohmic  potential 
drop. 

Srinivasan  et  al.  presented  experimentally  measured  cell 
voltage  data  obtained  at  368  K  with  H2/O2  at  4/5  atm,  as 
a  function  of  membrane  physico-chemical  characteristics 
(i.e.,  Nafion®  and  proprietary  Dow)  and  thickness  (i.e., 
llOOgeq-1  Nafion®  of  50, 100  and  175  pm  ‘dry’  thickness). 
The  membrane  ohmic  drop  was  estimated  from  the  cell  volt¬ 
age  using  a  proposed  cell  polarization  equation  and  afferent 
kinetic  parameters  [35].  Furthermore,  Meier  and  Eigenberger 
reported  the  membrane  ohmic  drop  as  a  function  of  current 
density  for  Nafion®  membranes  N112,  NE1135,  N115  and 
N1 17,  at  353  K  with  H2/02  at  2/2  bar  [31]. 

Figs.  2  and  3  show  the  dependence  of  7x5  on  1x2,  ex¬ 
tracted  from  the  data  reported  in  references  [35]  and  [31], 
respectively,  and  fitted  to  Eq.  (12)  extrapolated  to  7X2=0- 
The  hydraulic  pressure  drop  was  1  atm  in  Fig.  2  and  0. 1  atm 
was  assumed  for  the  conditions  relevant  to  Fig.  3.  The  cur¬ 
rent  density  ranges  were:  5000-14,000 Am-2  (Fig.  2)  and 
2000-12,500 Am-2  (Fig.  3),  respectively.  The  membrane 
physico-chemical  properties  are  summarized  in  Table  Al 
(Appendix  A).  The  literature  data  was  fitted  to  Eq.  (12)  writ¬ 
ten  as  7 r5  =  yo.m77!2- 

Figs.  2  and  3  show  that  the  ‘thin’  Nafion  membranes  ( ‘dry’ 
thickness  <  100  pm  such  as  N1 12)  and  the  Dow  membrane, 
were  characterized  by  a  linear  dependence  of  7x5  on  1x2-  How¬ 
ever,  for  Nafion  ‘dry’  thickness  of  127  pm  (N1 15,  Fig.  3)  and 


Fig.  2.  Correlations  between  dimensionless  numbers  7T5  and  712  for  selected 
Nafion  and  Dow  membranes.  Data  points  from  the  study  of  Srinivasan  et  al. 
[35].  Dry  membrane  thickness:  (•)  50  |xm  (Nafion),  (Q)  100  [xm  (Nafion), 
(▼)  125  |xm  (Dow),  (V)  175  |xm  (Nafion).  A p=  1  atm,  368  K. 

175-183  pm  (N117,  Figs.  2  and  3)  a  deviation  from  linearity 
was  obtained,  reflected  by  y2  of  1.3.  In  the  case  of  uniformly 
and  well-hydrated  membranes  a  linear  behaviour  is  expected 
for  the  membrane  potential  drop  as  a  function  of  current  den¬ 
sity  described  by  Ohm’s  law  [39].  The  latter  condition  and, 
therefore,  the  applicability  of  Ohm’s  law  is  likely  to  be  met 
at  low  current  densities  and  in  the  case  of  ‘thin’  membranes 
(e.g.,  up  to  about  100  pm  ‘dry’  thickness).  These  criteria  can 
be  better  expressed  using  the  dimensionless  number  1x2  (Eq. 
(7))  by  defining  a  linear-region  limiting  7T2,i  ratio.  Based  on 
Fig.  2,  7x2,1  =  0.035,  whilst  ;r2,i  =  0.02  in  the  case  of  Fig.  3. 
Therefore,  when  7x2  <  7x2,1  the  membrane  ohmic  drop  can  be 
estimated  with  A  <5  =  iZm/KQ. 

However,  for  7x2  >  txj,\  (i.e.  a  combination  of  high  cur¬ 
rent  density,  ‘thick’  membrane  and/or  low  fixed  charge-site 
concentration,  Eq.  (7))  the  ohmic  drop  behaviour  as  a  func- 


Fig.  3.  Correlations  between  dimensionless  numbers  715  and  712  for  selected 
Nafion  membranes.  Data  points  from  the  study  of  Meier  and  Eigenberger 
[31].  Nafion  membrane  types:  (•)  N1 12,  (O)  NE1 135,  (▼)  N1 15,  (V)  N1 17. 
A p  =  0.1  atm,  353  K. 
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Table  A1 

Membrane  properties 


Membrane 

type 

Typical  ‘dry’ 
thickness  (|xm) 

Basis  weight3 
(kgm^2) 

‘Wet’  thickness0 

Om) 

Acid  capacity 
(eq  kg~ 1 ) 

Fixed  charge-site  concentration  Cf, 
on  ‘wet’  membrane  basis  (molm-3)f 

N112 

51a 

0.10 

58 

0.89d 

1534 

NE1135 

OO 

VO 

0.19 

102 

idem 

1658 

N115 

127a 

0.25 

145 

idem 

1534 

N117 

183a 

0.36 

209 

idem 

1533 

Dow 

125— 130b 

NA 

148 

1.1 3b 

1918c 

a  Membrane  conditioned  at  296  K  and  50%  relative  humidity  (RH)  [36].  Note  the  Nafion  membranes  of  50,  100  and  175  |xm  dry  thickness  from  Ref.  [35] 
were  considered  to  have  similar  characteristics  as  N1 12,  NE1 135  and  N1 17,  respectively. 
b  Refs.  [35,37]- 

c  Estimated  from  the  ‘dry’  thickness  using  a  14%  thickness  increase  from  296  K  at  50%  RH  to  water  soaked  at  373  K  [36], 
d  From  Ref.  [36]. 

e  Estimated  assuming  a  25%  higher  fixed  charge-site  concentration  for  the  Dow  membrane  as  compared  to  N117  [5]. 

1  The  Cf  values  for  the  Nafion  membranes  are  close  to  the  average  of  1591  mol  m  1  given  in  Ref.  [38]  for  the  1 100  series  Nafion. 


tion  of  current  density  deviates  from  linearity.  This  is  further 
supported  by  a  separate  study  showing  the  increase  of  the 
experimentally  measured  in  situ  membrane  resistance  with 
current  density  for  membrane  ‘wet’  thickness  greater  than 
about  150  pun  (such  as  N117,  Appendix  A  and  Table  Al) 
[34].  Thus,  a  value  of  yz  >  1  in  Eq.  (13)  shows  that  the  mem¬ 
brane  resistance  increases  with  current  density,  possibly  due 
to  dehydration  of  some  parts  of  the  membrane  (most  likely 
on  the  anode  side)  during  the  operation  of  the  fuel  cell. 

Furthermore,  Eq.  (13)  estimates  fairly  well  the  membrane 
thickness  effect  on  the  ohmic  potential  drop.  In  the  experi¬ 
mental  study  of  Biichi  and  Scherer,  at  a  current  density  of 
5000 Am-2  for  N117  approximately  a  2.5  times  increase 
of  the  in  situ  membrane  resistance  was  projected  when  the 
‘wet’  thickness  of  N117  was  doubled  from  200  to  400  |xm 
[34].  According  to  Eq.  (13)  for  N1 17  with  yi  =  1.3,  doubling 
the  thickness  at  a  constant  current  density  implies  an  ohmic 
drop  increase  by  a  factor  of  about  2.5. 

It  is  interesting  to  note  at  this  point,  that  in  the  phenomeno¬ 
logical  modeling  of  the  Ballard  Mark  IV  fuel  cell  Amphlett 
et  al.  established  a  linear  dependence  for  the  internal  cell 
resistance  (due  mainly  to  membrane  N117  ohmic  loss)  as 
a  function  of  current  density  [19].  However,  the  same  au¬ 
thors  also  noted  that  quadratic  and  other  forms  of  non-linear 
effects  for  current  density  and  temperature  could  be  also  pos¬ 
sible.  Thus,  the  non-linear  dependence  of  membrane  ohmic 
loss  as  a  function  of  current  is  plausible  as  suggested  here 
for  N1 17. 

Regarding  the  exponent  yq  in  Eq.  (13),  which  reflects  the 
water  permeability  effect  on  the  membrane  ohmic  drop,  un¬ 
fortunately  neither  the  study  by  Srinivasan  et  al.  nor  the  one 
by  Meier  and  Eigenberger  do  not  present  enough  data  that 
could  be  used  to  estimate  independently.  However,  based 
on  the  work  of  Springer  et  al.  showing  that  the  high-frequency 
cell  resistance  under  well-humidified  conditions  is  indepen¬ 
dent  of  the  pressure  difference  A p  [40],  it  is  considered  that 
under  well-humidified  conditions  A<P  is  virtually  unaffected 
by  kpAp  (i.e.  the  streaming  potential  contribution  can  be  ne¬ 
glected).  The  well-humidified  membrane  condition  yields  V4 
equal  to  —  1 .  Dehydration  on  the  other  hand,  is  characterized 


by  y4<— 1,  since  an  increase  of  kpAp  would  decrease  the 
ohmic  potential  drop.  In  the  present  work  for  N117  it  was 
assumed  that  V4  =  —1.05,  allowing  for  approximately  a  10% 
decrease  of  A<£  for  a  10-fold  increase  of  kpAp. 

The  exponent  y(,  can  be  estimated  from  a  fit  of  7Z5  as  a 
function  of  jt(,  corresponding  to  different  open-circuit  condi¬ 
tions  of  humidity  and/or  temperature,  effecting  kq  (Eqs.  (10) 
and  (11)).  In  the  expression  of  7r(),  Eq.  (10),  the  fixed  surface 
charge  density  qf,  is  defined  per  geometric  cross-sectional 
area.  However,  it  is  more  relevant  to  express  qf  as  the  charge 
per  total  (or  effective)  membrane  surface  area.  The  latter  was 
calculated  using  the  equation  presented  by  van  der  Stegen  et 
al.  assuming  a  straight  and  parallel  cylindrical  pore  model  for 
the  membrane  [41]: 

(1  —  $m)kf|CfFrp 

*= - 2^ - ■  (14) 

where  rp  is  the  membrane  pore  radii  (m)  and  0m  the  water  up¬ 
take  of  the  membrane  per  volume  basis  (m3  m-3).  An  average 
pore  radius  of  2  nm  was  assumed  for  the  Nafion  membranes, 
based  on  the  study  by  Divisek  et  al.  indicating  that  the  lat¬ 
ter  average  pore  radius  corresponds  to  the  highest  fraction  of 
water  uptake  and  pore  volume  [42].  The  pore  size  of  the  Dow 
membrane  on  the  other  hand,  is  about  25%  smaller  compared 
to  Nafion  of  similar  thickness  [43],  therefore,  rp  =  1.5  nm. 

For  membranes  with  sulfonate  charge  sites,  9m  can  be 
estimated  by  [41]: 

9m  =  0.646  +  0.604 yw  +  0.1 12 y2  -  0.935EW 

+0.36EW2  -  0.441EWyw,  (15) 

where  EW  is  the  membrane  equivalent  weight  (kgeq-1, 
1.123  kg  eq-1  Nafion  and  0.885  kgeq-1  Dow,  see  Appendix 
A)  and  yw  is  the  water  activity  coefficient  (~0.9  to  0.975,  for 
1  =  9-14  [38]). 

From  Eqs.  (14)  and  (15)  for  the  conditions  explored 
in  Figs.  2  and  3  it  was  found  that  6*m  =  0.21  m3  m-3 
and  qf  =  0.55-0.6 Cm-2  for  the  Nafion  membranes  whilst, 
0m  =  O.45m3m-3  and  gf  =  0.17Cm-2  for  the  Dow  mem¬ 
brane,  respectively. 
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As  an  example  for  the  determination  of  yq  the  N1 17  mem¬ 
brane  was  employed.  The  jtq  ratio  between  the  conditions  of 
Figs.  3  and  2,  was  calculated  using  kq  =  13.2  Sm-1  (353  K 
and  X  =  16)  and  10.3  S  m-1  (368  K  and  X  =  11),  respectively. 
From  Eq.  (13),  using  the  775  ratio  between  Figs.  3  and  2  (i.e., 
equal  to  1 1.5  for  nx  =  const.)  and  the  previously  determined 
exponents,  yielded  y<,  =  —  1 .3  for  N1 17. 

Furthermore,  it  is  interesting  to  note  that  the  jte  values 
were  between  0.010  and  0.016  for  the  Nafion  1100  series  and 
0.23  for  the  Dow  membrane,  respectively.  These  ratios  indi¬ 
cate  that  the  actual  membrane  conductivity  for  Nafion  is  up  to 
two  orders  of  magnitude  smaller  than  the  surface-convective 
conductivity  supporting,  therefore,  the  hypothesis  that  the 
membrane  proton  conductivity  is  not  predominantly  surface 
mechanism-based  [44].  However,  the  approximately  twenty 
times  higher  rtf,  ratio  for  the  Dow  membrane  compared  to  the 
Nafion,  suggests  that  in  the  case  of  the  former  the  surface- 
convective  mechanism  plays  a  more  significant  role. 

Based  on  the  proposed  dimensionless  groups  and  mono¬ 
mial  correlating  function,  for  1x2  >  Tt2}  the  ohmic  drop  for 
the  N117  membrane  under  the  conditions  presented  in 
Figs.  2  and  3  can  be  expressed  as 


\zi\FCf(k0A<P)\ 


kvAp 


-) 


=  5  x  10" 


—3  (  t  Zm-Sji+P 


nF  |zf I  Cfp. 


1.3 


(kpAp)p\ 
F2  ) 


-1.05 


AT0/X 

,  f  , 


-1.3 


(16) 


Note  that  as  discussed  previously,  in  Eq.  (16)  the  specific 
conductivity  at  open-circuit  conditions  is  used,  which  is  more 
readily  accessible  experimentally  and/or  by  calculation  as 
compared  to  the  in  situ  specific  conductivity. 

The  coefficient  yo.m  in  Eq.  (16),  can  be  interpreted  as  the 
value  of  7T5  in  the  particular  case  when  772  =  774  =  776  =  1 .  Since 
typical  values  of  774  and  Jt(,  are  about  10-4  to  10  ’  and  1 0  2 
to  10“ 1 ,  respectively,  the  special  case  of  774  =  ne  =  1  indicates 
a  combination  of  high  hydraulic  water  permeability  and  ionic 
conductivity,  therefore,  it  is  expected  that  yo,m<  1  (see  also 
Eq.  (9)). 


4.2.  Dimensionless  numbers  and  monomial  correlating 
equations  for  the  catalyst  layer 

Table  2  shows  the  dimensionless  numbers  describing  the 
catalyst  layer.  The  exponential  term  containing  the  electrode 
potential  (0S  —  0)  was  assigned  to  the  dimensionless  num¬ 
ber  7ro .  Equivalent  forms  for  ttq  can  be  written  in  terms  of 
overpotential  q  =  0S  —  0  —  £°  (as  defined  in  Ref.  [39])  and 
Tafel  slope  for  O2  reduction  b  =  23RT/aF: 

aF  \ 

~rRTT]) 


77q  =  exp (~f(0s  -0  -  ED)  =  exp 


=  exp 


where  commonly  the  transfer  coefficient  a  =  0.6  —  1  based 
on  experimentally  measured  Tafel  slopes  for  O2  reduction 
in  acid  [39,40].  For  negligible  ohmic  potential  drop  in  both 
the  electronic  (solid-matrix)  and  ionic  (polymer  electrolyte) 
conductive  phases  of  the  catalyst  layer,  the  overpotential  in 
Eq.  (17)  is  constant  across  the  catalyst  layer  thickness  [4,32]. 

The  dimensionless  number  779  (Table  2)  compares  the  in¬ 
trinsic  kinetic  rate  of  the  electrode  reaction  with  the  effective 
rate  of  O2  diffusion  in  the  catalyst  layer  (Eq.  (18)).  Therefore, 
779  could  be  considered  the  electrode  kinetic  equivalent  of  the 
second  Damkohler  number  Dan  (for  a  review  see  [22]).  Alow 
Da\\  in  the  present  case  means  that  the  effective  O2  diffusion 
rate  is  high  across  the  catalyst  layer  thickness  compared  to 
the  intrinsic  electrode  kinetic  rate  defined  at  £°  (note  the  k° 
inEq.  (18)). 

Furthermore,  it  is  easily  recognized  that  the  square  root  of 
779  is  similar  to  the  expression  of  the  Thiele  modulus  4>  for 
first-order  kinetics  [45]: 

(. ack°)Z ;  7 

779  =  Dan  =  n  c  =  02,  (18) 

D02.C 

where  ac  is  the  specific  surface  area  of  the  catalyst  layer, 
Do2,c  the  effective  O2  diffusion  coefficient  in  the  catalyst 
layer,  k°  the  standard  heterogeneous  rate  constant  and  Zc  the 
catalyst  layer  thickness. 

The  Thiele  modulus  as  defined  by  Eq.  (18)  is  related  to 
the  standard  heterogeneous  rate  constant.  However,  it  is  im¬ 
portant  to  consider  also  the  electrode  potential  dependency 
by  employing  the  product  779770.  Thus,  from  Eqs.  (17)  and 
(18)  an  overpotential-dependent  Thiele  modulus  <j>E  can  be 
defined  as 


JTgJTQ  = 


<Pi 


(ack°)Z2c 

D02.C 


exp 


=  Da  1 1  exp 


(19) 


The  electrode  potential-dependent  Thiele  modulus  (pE  = 
779770,  reflects  the  catalyst  layer  utilization  efficiency  and 
penetration  depth.  It  must  be  noted  that  Perry,  Newman  and 
Cairns  in  the  analysis  of  liquid-electrolyte  fuel  cells,  intro¬ 
duced  a  Thiele  modulus  conceptually  similar  to  <Pe  but  ex¬ 
pressed  in  terms  of  the  porous  catalyst  agglomerate  radius 
instead  of  Zc  [46].  The  variant  given  in  the  present  work, 
on  the  other  hand,  is  based  on  the  pseudohomogeneous  film 
model  of  the  catalyst  layer. 

From  Eq.  (18)  for  Zc=  const.,  the  higher  (p  or  Dan  the 
smaller  the  penetration  depth  of  the  O2  reduction  front  inside 
the  catalyst  layer.  Furthermore,  from  Eq.  (19),  for  the  same 
Dan  the  more  negative  the  cathode  overpotential,  the  higher 
<pE  and  consequently,  the  lower  the  penetration  depth  of  the 
reaction  front.  For  efficient  utilization  of  the  Pt-based  cata¬ 
lysts  the  physical  thickness  of  the  catalyst  layer  Zc  should  be 
approximately  equal  to  the  penetration  depth  up  given  by  the 
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following  equations: 


_  Zc  _  f  %,c  1/2 

Vp  cp  v  ack°  ) 


Vp,E  = 


Zc 

0E 


{  %,c 

V  flc  k° 


(20) 

(21) 


Eqs.  (20)  and  (21)  define  the  penetration  depth  in  terms  of 
electrode  kinetics  coupled  with  O2  mass  transport  based  on 
the  pseudohomogeneous  film  model.  Newman  and  Tiede- 
mann  defined  a  porous  electrode  penetration  depth  for  the 
case  of  electrode  kinetics  in  conjunction  with  ohmic  limita¬ 
tion  [24].  Their  ratio  was  also  inversely  proportional  to  the 
square  root  of  the  specific  surface  area-exchange  current  den¬ 
sity  product  acio- 

The  rest  of  the  dimensionless  numbers  in  Table  2,  i.e.  from 
7T 10  to  jt  1 7,  can  be  related  to  dimensionless  numbers  that 
have  already  been  introduced.  Thus,  7Tio  (Eq.  (22))  could  be 
considered  the  Nusselt  number  for  O2  transport  in  the  cat¬ 
alyst  layer  (see  also  Eq.  (8)  and  Ref.  [32]),  assuming  diffu¬ 
sion  across  the  entire  catalyst  layer  thickness  and  fast  elec¬ 
trode  reaction  concentrated  at  catalyst  layer/bulk  membrane- 
electrolyte  interface  (Fig.  1)  with  ACo2,c  =  Co,  1  —  Co,2> 
where  Cq2,2  =  0  and  ACq2  =  Cq2,i  =  Cq2 : 


trio  = 


/  Zcv< ), 

nFDo2XCo2 


i 

!l,o2,c 


(22) 


The  concentration  C o2,  corresponding  to  the  catalyst/gas- 
diffusion  layer  interface,  can  be  calculated  from  Henry’s  law: 
C’0,1  =  A'O),  1  Pd/ H  with  pi  the  inlet  air  pressure  in  the  dif¬ 
fuser,  xq,  |  the  O2  mole  fraction  at  the  catalyst/gas-diffusion 
layer  interface  and  //is  Henry’s  law  constant  [4], 

The  numbers  t*\h  and  7115  on  the  other  hand,  express  the 
ohmic  drop  in  the  membrane -electrolyte  A<P  and  solid  phases 
A<Z>S  of  the  catalyst  layer,  respectively,  similarly  to  7*5  but 
taking  into  account  the  effective  values  in  the  porous  media 
of  the  respective  physico-chemical  constants,  k$  c  and  kpx. 

The  dimensionless  numbers  7ri6  and  tt  1 7  are  similar  to 
it(,  (Table  1).  Thus,  7Ti6  and  jr  1 7  compare  the  effective  ionic 
kc  and  electronic  erc  conductivities,  respectively,  with  the 
surface-convective  proton  conductivity  of  the  membrane- 
electrolyte: 


ionic  transport  (i.e.  ionic  conductivity  effect)  and  (IV)  elec¬ 
trode  kinetics,  O2  and  ionic  transport  control.  It  is  proposed 
that  the  catalyst  layer  polarization  curves  corresponding  to 
the  various  operating  regimes  can  be  described  in  a  dimen¬ 
sionless  form  using  a  monomial  correlating  function  and  4>p 
as  the  dependent  variable: 

06  =  (25) 


or  explicitly. 


(ack°)Z2c  (  2.3  \ 

£>02,c  6XP  V  b  V 

{  i Zcsq2  \ no  (  ^cm42  \ m 

'  °  c  \nFD02XCo2)  \^{n\zi\FCf)2 Z2  J 


kfl  Cf 
Co2 


y  12 


(26) 


To  reduce  the  number  of  unknown  exponents,  j*k,  and  t*\2 
could  be  combined  in  one  term,  7ri6 ,c,  considering  y  12  =  2_v’u,, 
which  leads  to  the  following  equations: 


</>£ 


\jr\  7T^10  7T'^i6 

yO,cK  10  71 16, C 


(27) 


and 


(ack°)Z2c  f  2.3  \ 

Dq2,c  6XP  \  b  7 

/  izcs02  yiol  kc^2q2  \m 
0c\nFDo2XCo2)  \(nFCo2)2Z2J 


(28) 


For  regime  I:  intrinsic  electrode  kinetic  control,  the  ef¬ 
fective  O2  diffusivity  and  ionic  (i.e.  proton)  conductivity  do 
not  affect  the  performance,  therefore,  in  Eqs.  (27)  and  (28) 
simply  yi6  =  0  while  yio  =  1,  with  yo,c  =  1  one  obtains: 

4>l  =  nio  (29) 

and 

n  n  /  2.3  \ 

1  =  — FZcack°Co2  exp  -  —  rj  .  (30) 

^o2  V  b  ) 


7*16 


KcFSq2 

(n\z{\FC{)2  Z2 


K  cM 


and 


7*11 


(n\zf\FCf)2Z2 


acp 


(23) 


(24) 


Regarding  the  polarization  behaviour  of  the  catalyst  layer, 
the  following  operating  regimes  can  be  identified  [46]:  (I)  in¬ 
trinsic  electrode  kinetic  control,  (II)  mixed;  electrode  kinet¬ 
ics  and  O2  mass  transport  control,  (III)  kinetics  coupled  with 


Moreover,  y  10  was  also  determined  from  a  h  t  of  i]  as  a  func¬ 
tion  of  In  i,  with  the  slope  equal  to  —yio/V2.3  from  Eq.  (29) 
written  as  =  Tr-}),11 .  The  experimental  cathode  polarization 
curve  at  5  atm  neat  O2  pressure  reported  by  Springer  et  al. 
(Fig.  15,  curve  a,  in  Ref.  [40])  was  employed.  For  the  current 
density  range  of  740-8000  Am-2,  using  6  =  0.085  V dec-1 
[40],  it  was  found  that  yio  =  1,  confirming,  therefore,  the  de¬ 
pendence  expressed  by  Eqs.  (29)  and  (30). 

In  the  case  of  mixed  control  according  to  regime  II:  elec¬ 
trode  kinetics  coupled  with  O2  transport,  using  yi6  =  0  and 
Vo,c  =  1  -  the  general  equation  is 

06  =  ^io0 


(31) 
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or 

i  =  —  FCO2Zc({2/m)-l\ack0)l/m 
502 

x  (£>O2,c)(1“(1/),10,)exp  •  (32) 

The  general  equation  for  regime  II  (i.e.,  Eq.  (32))  shows 
the  characteristic  multiple  Tafel  slope  behaviour  of  the  porous 
electrode: 

b*  =  ywb.  (33) 

Eqs.  (31)  and  (32)  show  the  power  of  dimensional  analysis 
to  provide  in  a  mathematically  simple  way  a  general  analyti¬ 
cal  solution  to  a  complex  problem.  Employing  the  experimen¬ 
tal  cathode  polarization  curve  obtained  by  Springer  and  co¬ 
workers  at  5  atm  air  pressure  (Fig.  15,  curve  b,  in  Ref.  [40]), 
y  i  o  was  determined  from  a  fit  to  Eq.  (31),  as  described  previ¬ 
ously.  Using  the  current  density  range  of  3000-9100  Am-2, 
yielded  yio  =  2.  Thus,  according  to  Eq.  (33)  a  double  Tafel 
slope  was  obtained  for  the  catalyst  layer  operating  on  5  atm 
air,  i.e.  b  =  2b.  The  same  Tafel  slope  was  also  obtained  us¬ 
ing  the  data  for  2  atm  O2/N2  mixture  with  13.5  vol.%  O2  for 
current  densities  up  to  5000  A  m-2  [40].  Eq.  (32)  in  the  case 
of  yio  =  2  becomes 

n  n  1/2  (  2.3  \ 

i =  — FCo2(ack0Do,x )  exP  '  (34) 

so2  V  2  b  J 

Using  an  asymptotic  approach  [47],  Perry  et  al.  derived 
equations  for  the  catalyst  layer  under  kinetic  and  mixed,  ki¬ 
netic  and  O2  mass  transport  control,  respectively,  similar  to 
Eqs.  (30)  and  (34)  (i.e.  particular  case  of  vio  =  2)  [46].  Com¬ 
paring  Eqs.  (30)  and  (34)  the  double  Tafel  slope  is  indicative 
of  mass  transfer  effects  limiting  the  current  density,  while  for 
both  cases  the  oxygen  concentration  dependence  is  of  first 
order. 

In  addition  to  kinetics  and  O2  mass  transfer,  the  ionic  con¬ 
ductivity  could  affect  the  performance  of  the  catalyst  layer. 
Experimental  results  [40]  indicated  that  current  ratio  mea¬ 
sured  at  0.7  V  between  the  runs  with  5  atm  air  and  2  atm 
O2/N2  mixture  with  13.5%  O2,  was  around  2  instead  of  about 
5,  what  would  have  been  expected  based  on  the  first-order  O2 
partial  pressure  dependence  (Eq.  (34)).  This  limitation  was 
attributed  to  ionic  conductivity  effects  [40].  The  latter  could 
be  accommodated  by  using  the  ionic  conductivity  term  in 
Eqs.  (27)  and  (28)  and  determining  yi6  in  addition  to  the  O2 
mass  transfer  effect  characterized  by  yio  =  2  as  shown  before. 
Thus,  based  on  Eq.  (28)  with  yio  =  2,  the  current  density  is 
proportional  to  the  O2  concentration  according  to: 

*~Cgf+1.  (35) 

Substituting  in  Eq.  (35)  the  cathode  current  densities  mea¬ 
sured  atO. 7  V (Fig.  15, Ref.  [40]), i.e.,  ll,500Am-2  for5  atm 
air  and  6000  A  m-2  for  2  atm  O2/N2  mixture  with  13.5%  O2, 


Fig.  4.  Catalyst  layer:  overpotential  as  a  function  of  current  density 
for  various  oxygen  reduction  rate-limiting  scenarios.  Summary  of  condi¬ 
tions:  Zc  =  10  (Jim,  kc  =  0.25  S  m-1,  b  =  0.085  V  dec-1,  ac  =  2  x  107  m2  m-3, 
D0,  c  =  2  x  10-8  m2  s-1,  T=353K. 

a  half-order  O2  concentration  dependence  of  the  current  den¬ 
sity  is  obtained.  Hence  yi6  =  —0.5  and  Eq.  (27)  becomes 

4  =  *W,c-  (36) 

From  Eq.  (36)  the  polarization  equation  in  the  case  of 
combined  ionic  and  O2  mass  transport  control  coupled  with 
electrode  kinetics  is 

i  =  ^ 5-Fack°Do2,c(Y1(KcH)1/2  exp  ^-^77^  •  (37) 

To  the  knowledge  of  the  author  equations  (36)  and  (37)  are 
the  first  analytical  equations  describing  the  proton  transport 
effect  combined  with  kinetics  and  O2  transport  (regime  IV). 
The  current  density  shows  a  half-order  dependence  on  the 
O2  concentration  as  observed  experimentally  [40].  Moreover, 
Eq.  (37)  gives  a  1/4-order  dependence  of  i  on  kc,  which  pre¬ 
dicts  fairly  well  the  current  density  ratios  reported  by  Springer 
et  al.,  corresponding  to  0.7  V  cathode  potential  and  effec¬ 
tive  ionic  conductivities  of  the  catalyst  layer,  of  0.5,  0.2  and 
0.08  S  m”1  (Fig.  7  in  Ref.  [40]). 

Eq.  (37)  shows  that  in  the  case  of  both  O2  and  ionic  trans¬ 
port  limitation,  the  current  density  is  inversely  proportional 
with  the  square  root  of  catalyst  layer  thickness  for  a  constant 
catalyst  specific  surface  area  ac.  Therefore,  a  thin  catalyst 
layer  is  required  to  reduce  the  impact  of  the  limiting  trans¬ 
port  steps.  Finally,  Eq.  (37)  with  ;/  constant  across  the  catalyst 
layer  thickness  implies  that  the  ionic  conductivity  effect  re¬ 
flects  not  an  ohmic  control  but  a  proton  reactant  (Eq.  (2)) 
diffusion  limitation. 

Fig.  4  compares  the  effects  of  the  various  lim¬ 
iting  factors  on  the  catalyst  layer  polarization  be¬ 
haviour  for  5  atm  air  at  353  K  according  to  Eqs.  (30), 
(34)  and  (37),  respectively,  assuming  no  gas  diffusion¬ 
backing  layer  limitation.  The  physico-chemical  parame¬ 
ters  were:  Zc=10|xm,  kc  =  0.25  S  m“ 1 ,  (7  =  0.085  V  dec-1 
[40],  ac  =  2x  107  m2  m-3,  Dq2,c  =  2  x  10“8  m2  s”1  [39] 
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Fig.  5.  Catalyst  layer:  Thiele  modulus  for  electrode  kinetics  0e  as  a  function 
of  trio-  Conditions  identical  to  Fig.  4. 


and  Co  =  4. 96  mol  m-3  corresponding  to  5  atm  air  and 
Henry’s  constant  of  2.04  x  104  Pam3  mol-1  calculated  from 
the  formulae  given  by  Bernardi  and  Verbrugge  [4]  for  the  O2 
solubility  in  Nafion.  The  standard  heterogeneous  rate  con¬ 
stant  A:0  was  7.9  x  10-11  ms-1,  calculated  from  Eq.  (3)  using 
io=  1.5  x  10-4  Am-2  and  Co2,ref  =  4.9 mol m-3  for  lbar 
O2  pressure  [39]. 

Additionally,  the  combined  influence  of  electrode  kinetics 
and  ohmic  drop  (i.e.  assuming  fast  O2  transport)  was  plotted 
using  the  equation  derived  by  Perry  et  al.  based  on  asymptotic 
analysis  [46]: 


/ 2«o, Fack^Co^Kcb  f  2.3  \ 

' =  V  23  exp  \2b")-  (38) 

Fig.  5  on  the  other  hand,  shows  the  same  effects  as  Fig.  4 
but  expressed  in  dimensionless  form  as  the  dependence  of  the 
overpotential-dependent  Thiele  modulus  <pE  as  a  function  of 
the  O2  Nusselt  number  for  the  catalyst  layer,  tz\$. 

As  shown  in  Fig.  4,  superficial  current  densities  up  to 
23,000 Am-2  could  be  achieved  in  the  case  of  intrinsic 
electrode  kinetic  control  (Eq.  (26))  with  cathode  overpo¬ 
tentials  —  ?7<0.5V.  Both  O2  and  proton  transport  limita¬ 
tions,  separately  and  in  combination,  have  the  role  of  in¬ 
creasing  the  overpotential  required  for  a  given  superficial 
current  density.  The  combined  limiting  effect  of  O2  and 
ionic  transport  increases  dramatically  the  cathode  overpo¬ 
tential,  e.g.  a  current  density  of  about  3000 Am-2  was 
obtained  at  an  overpotential  of  —0.7  V  (Fig.  4).  In  di¬ 
mensionless  form.  Fig.  5  shows  that  at  constant  jrio  the 
overpotential-dependent  Thiele  modulus  was  much  higher 
when  both  ionic  and  O2  transport  limited  the  electrode  per¬ 
formance. 

The  penetration  depth  vPie  according  to  Eq.  (21)  at 
7Tio  =  0.6  (i.e.  2300  Am-2)  was  about  15  pm  for  O2  trans¬ 
port  limitation  and  0.4  pm  for  both  ionic  and  O2  transport 
control.  This  result  shows  clearly  the  inefficient  utilization 


of  the  catalyst  layer  thickness  in  the  latter  scenario  (physical 
thickness  of  10  pm). 

Based  on  Fig.  5  the  parameter  r  =  could  serve  as  an 
indicator  for  the  different  operating  regimes: 

•  Regime  I:  intrinsic  electrode  kinetics  r  <  1, 

•  Regime  II:  electrode  kinetic  and  O2  mass  transport  r  =  1, 

•  Regime  III:  electrode  kinetic  and  ionic  transport  r  >  1,  and 

•  Regime  IV:  electrode  kinetics  with  O2  and  ionic  transport 
r  »  1. 

The  ohmic  potential  drop  in  both  the  electrolyte  and  elec¬ 
tronic  conductor  phases  of  the  catalyst  layer  can  be  expressed 
with  the  dimensionless  ratios  given  in  Table  2.  For  the  elec¬ 
trolyte  phase  the  correlating  function  is  similar  to  the  one  de¬ 
scribed  for  the  bulk  electrolyte-membrane  (see  Section  4.1) 
but  in  the  expression  of  the  dimensionless  ratios  the  effective 
values  of  the  physico-chemical  parameters  are  used.  In  the 
case  of  the  electronic  conductor  solid  phase.  Ohm’s  law  in  di¬ 
mensionless  form  is  written  as:  7715  =  jtioJtfjtflTtfjtf). 
Furthermore,  depending  on  the  effective  conductivities  of  the 
two  phases  crc  and  kc,  the  ohmic  potential  drop  for  the  elec¬ 
trolyte  and/or  solid  phase  can  be  included  in  the  total  po¬ 
larization  behaviour  of  the  catalyst  layer  [24,46].  One  way 
to  evaluate  the  importance  of  the  ohmic  drop  is  to  compare 
the  charge-transfer  and  ohmic  resistances  by  calculating  the 
Wagner  numbers  for  the  catalyst  layer,  Wa  =  Kc(Sq/Si)/Zc  and 
Was  =  ac(Sri/Si)/Zc,  using  the  polarization  expressions  given 
by  Eqs.  (30),  (34)  and  (37).  For  Wa,  Was  <  1  the  ohmic  resis¬ 
tance  dominates  and  it  must  be  accounted  for  in  the  polariza¬ 
tion  equation. 


4.3.  Dimensionless  numbers  and  monomial  correlating 
equations  for  the  gas-diffusion  (backing)  layer 


Table  3  shows  the  dimensionless  numbers  derived  for 
the  gas-diffusion  layer  using  the  Quraishi-Fahidy  technique. 
One  goal  for  the  gas-diffusion  layer  analysis  is  to  evalu¬ 
ate  the  effect  of  mass  transport  on  the  O2  mole  fraction 
at  the  catalyst/gas  diffusion-backing  layer  interface,  xo2,i 
(see  Fig.  1).  Therefore,  in  Table  3  the  O2  mole  fraction 
difference  across  the  diffuser  has  been  employed,  Axo,  = 
xo,,0  ~  To  1 ,  where  xo,o 's  the  inlet  O2  mole  fraction. 

Relevant  mass  and  momentum  transport  related  dimen¬ 
sionless  numbers  characteristic  for  the  gas-diffusion  layer  are 
(Table  3):  Schmidt  numbers  for  the  effective  binary  mixture 
diffusivities  of  O2-N2  7Ti8,  and  water-air  7ri9,  respectively, 
the  water  Reynolds  number  in  the  diffuser  nn,  and  the  Nus¬ 
selt  number  for  O2  mass  transfer  in  the  gas  diffusion  layer 
7T23  or  jt23,w  (see  also  Eq.  (8)): 


7T23 


isSQ2ZdRT 
n  TDn2— o2,d  Axq2  pd 


Nu  E,d, 


(39) 


^23, w 


isso2  Z&RT 
nFDy,— air,  dAxo2  Pd 


/Vt<E,d,w> 


(40) 
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where  pd  is  the  inlet  air  pressure  at  the  cathode  diffuser, 
D\2_02:d  the  effective  binary  O2-N2  diffusivity  (a  function 
of  wet-proofed  pore  volume  fraction),  Z)w_air,d  the  effective 
water-air  binary  diffusivity  (a  function  of  the  total  pore  vol¬ 
ume  fraction  in  the  absence  of  wet-proofing),  is  the  superficial 
current  density  in  the  diffuser  solid  phase  (equal  to  the  oper¬ 
ating  current  density  it)  and  Zd  the  diffusion  layer  thickness. 

Eqs.  (39)  and  (40)  are  based  on  the  pseudohomogeneous 
film  model  considering  that  the  electrochemical  reaction  is 
concentrated  at  the  catalyst/gas -diffusion  layer  interface  (po¬ 
sition  1,  Fig.  1).  Moreover,  Eqs.  (39)  and  (40)  distinguish 
between  two  scenarios  as  suggested  by  the  respective  diffu- 
sivities.  Eq.  (39)  considers  a  separate  network  of  hydrophobic 
and  hydrophilic  pores  in  the  diffuser  (realized  in  practice  by 
addition  of  a  wet-proofing  agent  such  as  PTFE)  and  con¬ 
sequently,  air  diffuses  exclusively  through  the  hydrophobic 
pore  network.  While  Eq.  (40),  by  employing  Z?w_air,c 1  inher¬ 
ently  implies  that  ah'  and  water  share  the  same  network  of 
pores. 

The  mass  transfer  limiting  condition  in  Eqs.  (39)  and  (40) 
is  given  by  x'0,,1  =  0  (i.e.,  the  O2  is  completely  consumed 
by  fast  electrode  reaction  at  the  catalyst/gas-diffusion  layer 
interface).  This  leads  to 


tf23.L 


isso2ZdRT 

nFDN2-o2zdxo2loPd 


is 

'L,02,d 


and 


^23,L,w 


issp2ZdRT 

11 F /4W— air,  dV()2.0  Pd 


*L,02,d,w 


(41) 


(42) 


where  tL,o2,d>  and  T.Oixtw  are  the  gas  diffuser  mass  transfer 
limited  superficial  current  densities  without  and  with  inter¬ 
ference  from  a  co-diffusing  water  phase,  respectively. 

A  monomial  correlating  equation  for  the  case  of  shared 
pore  network  can  be  written  as 


^23,  w  =  yo.d^if  i-e-’  N«E,d,w  =  yo.d^22.  (43) 


A  high  Rc\v  in  the  diffuser  describing  water  flow  from 
the  catalyst  layer  into  the  gas  diffusion  layer  (i.e.,  flooding) 
impedes  the  O2  gas  transport  and  reduces  Nue4,w  in  other 
words,  Axq2  =  A'0,,0  —  -A),.  1  increases.  Therefore,  V22  < 9 
in  the  case  of  backing-diffusion  layer  flooding. 

Substituting  the  corresponding  dimensionless  ratios  in  Eq. 

(43)  gives  xo2,i  in  the  scenario  of  shared  pore  network  for 
both  air  and  water: 

*o2,i  _  j _ 1_  /  isso2ZdRT  \  / Up  dpZd \  ~-'-2 

-*O2,0  >’0,d  \  n  TOw— air,d-^O2,0Pd  /  \  H  ) 

=  1-—  7r23,L,w/?d22.  (44) 

yo,d 

As  an  example,  using  the  superficial  water  velocities  reported 
by  Bernardi  and  Verbrugge  at  4000  and  6000  Am-2  (Fig.  9, 
Ref.  [5])  and  assuming  Axo2  =  0.05xo2,o  for  the  former  and 
Axo,  =  0.  lxo2,o  for  the  latter  superficial  current  density,  Eq. 

(44)  yields  V22  =  —0.45  and  yo,d  =  10. 


In  the  case  of  separate  hydrophobic  and  hydrophilic  pore 
networks  in  the  diffuser  maintained  throughout  the  operation 
of  the  fuel  cell,  Eq.  (44)  reduces  to 


xo2,i 

xo2,o 


1 

=  1 - 7T23.L  =  1  “ 

>’0,d 


1  is 

VO.d  'L.Oo.d 


(45) 


Substituting  the  mole  fraction  ratios  for  1000,  5000  and 
8000  Am”2  from  Fig.  7  of  Ref.  [4],  yields  y0,d=  1.  Thus, 
Eq.  (45)  is  consistent  with  the  familiar  form  of  the  surface 
concentration  obtained  using  Fick’s  law. 

The  ohmic  potential  drop  in  the  solid  electronic  conductor 
phase  can  be  expressed  in  a  dimensionless  form  as 


*2A  =  yQ,dn$n$n$Jt)gi,  (46) 

where  d  =  1,  V23=  1  and  Vis  =y22  =>’25  =  - 1 ,  yields  the 
expression  for  Ohm’s  law. 


4.4.  Applications 

Sensitivity  analysis  is  performed  on  the  developed  corre¬ 
lating  equations  with  respect  to  key  variables  and  exponents 
such  as  yio  in  Eq.  (32)  controlling  the  multiple-Tafel  slope 
value  of  the  mixed:  kinetic  and  O2  mass  transfer  regime,  V22 
in  Eq.  (43)  indicative  of  diffusion-backing  layer  flooding  and 
the  combined  influence  of  catalyst  layer  thickness  and  ef¬ 
fective  conductivity  on  regime  IV  (Eq.  (37)).  Furthermore, 
the  dimensional  analysis-based  model  predictions  including 
the  membrane  ohmic  drop  described  by  Eq.  (13),  are  com¬ 
pared  with  experimental  data  independent  of  the  those  used 
to  develop  these  correlations. 

Fig.  6  shows  the  effect  of  yio  on  the  cathode  potential 
Eoc  +  )]  for  air  pressure  of  5  atm,  at  353  K.  It  was  assumed 
that  the  open-circuit  value  Eoc  =  0.98  V  due  to  the  establish¬ 
ment  of  a  mixed  potential  on  the  cathode  surface.  A  separate 
hydrophilic-hydrophobic  pore  network  was  considered  in  the 
diffuser  expressed  by  Eq.  (45)  (i.e.  no  gas  diffuser  flooding 
limitation).  The  catalyst  layer  thickness  was  25  |xm  and  the 
specific  surface  area  of  the  catalyst  2  x  109  m2  m-3.  The  rest 


Fig.  6.  Effect  of  parameter^  10  on  the  catalyst  layer  polarization  under  mixed, 
electrode  kinetic  and  O2  mass  transfer  control.  5  atm  air  pressure,  353  K. 
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Fig.  7.  Effect  of  diffuser  flooding  on  the  catalyst  layer  polarization  under 
mixed,  electrode  kinetic  and  O2  mass  transfer  control  with  yio  =  2.  5  atm  air 
pressure,  353  K. 

of  the  physico-chemical  data  was  given  in  Section  4.2.  The 
dimensionless  number  jtio  could  serve  as  an  indicator  of  in¬ 
trinsic  kinetic  or  mixed  (kinetic  and  Ot  transport)  control.  Un¬ 
der  the  conditions  of  Fig.  6,  Eqs.  (30)  for  kinetic  and  (32)  for 
mixed  control  (regime  II),  respectively,  intersect  at  7Tio  =  0.4. 
Regime  II  dominates  at  7Tio>0.4  with  y  i  o  (multiple-Tafel 
slope  value)  affecting  the  polarization  behaviour,  whilst  for 
jtio  <  0.4  the  polarization  behaviour  is  described  by  regime  I 
(i.e.  intrinsic  electrode  kinetic  control).  Fig.  6  shows  that  the 
transition  between  the  two  regimes  occurred  at  1500  Am-2. 

In  regime  II  the  O2  penetration  depth  becomes  smaller  than 
the  physical  thickness  of  the  catalyst  layer.  In  other  words,  the 
y  10  effect  is  a  reflection  of  the  different  oxygen  permeabilities 
of  the  catalyst  layer. 

Fig.  7  shows  the  effect  of  diffuser  flooding  on  the  polariza¬ 
tion  behaviour  using  Eq.  (44)  to  calculate  the  O2  concentra¬ 
tion  at  the  diffuser/catalyst  layer  interface.  Three  yji  values 
were  compared,  —0.6,  —0.65  and  —0.7,  respectively,  whilst 
yo,d  =  10  (Section  4.3)  andyio  =  2.  The  diffuser  thickness  was 
260  p,m.  The  rest  of  the  physico-chemical  parameters  were 
the  same  as  presented  previously.  Fig.  7  reveals  that  the  dif¬ 
fuser  flooding  can  be  succinctly  described  by  Eqs.  (43)  and 
(44)  coupled  with  Eq.  (32)  for  regime  II.  Thus,  the  proposed 
monomial  correlations  eliminate  the  need  for  the  exponential 
term  used  previously  in  the  literature  to  model  empirically  the 
O2  mass  transfer  (Eq.  (1))  [18].  The  more  negative  V22  the 
more  severe  the  diffuser-backing  layer  flooding,  compromis¬ 
ing  the  cathode  performance.  Under  the  conditions  explored 
in  Fig.  7  at  5000  A  m-2  the  cathode  potential  in  the  absence 
of  flooding  (i.e.  separate  hydrophilic  and  hydrophobic  pore 
network)  was  0.67  V,  while  0.62  and  0.47  V  were  obtained  in 
the  case  of  flooding  characterized  by  >'22  =  —0.65  and  —0.7, 
respectively  (Fig.  7). 

The  proton  transport  influence  on  the  cathode  polariza¬ 
tion  behaviour  can  be  detected  by  the  combined  effect  of 
the  catalyst  layer  thickness  Zc  and  specific  protonic  conduc¬ 
tivity  kc.  Assuming  vio  =  2,  Eq.  (37)  shows  that  in  the  case 


Fig.  8.  Comparison  between  regime  II  (kinetic  andOa  mass  transfer  control) 
and  regime  IV  (kinetic,  O2  and  ionic  mass  transfer  limitation)  for  catalyst 
layer  polarization.  5  atm  air  pressure,  353  K,  catalyst  layer  thickness:  5  and 
25  |xm,  effective  ionic  conductivity  of  the  catalyst  layer  0.25  and  2.5  S  m  1 . 

of  both  ionic  and  O2  transport  limitation  (regime  IV)  for  a 
given  overpotential  the  current  density  is  proportional  to  the 
square  root  of  the  catalyst  layer  thickness,  whilst  for  intrinsic 
O2  mass  transport  control  (Eq.  (34),  regime  II),  the  polariza¬ 
tion  behaviour  is  independent  of  thickness.  Fig.  8  compares 
the  two  polarization  regimes  II  and  IV,  in  the  case  of  5  and 
25  |xm  catalyst  layer  thickness  and  0.25-2.5  S  m-1  protonic 
conductivity,  respectively.  It  was  assumed  that  the  cathode 
is  not  limited  by  diffuser  flooding  and  the  specific  catalyst 
surface  area  is  the  same  in  all  cases.  Under  the  conditions  of 
Fig.  8,  the  combined  limitation  of  proton  and  O2  transport 
in  the  catalyst  layer  decreases  the  cathode  potential  between 
about  0.2  and  0.6  V  at  current  densities  above  1000 Am-2. 
Based  on  Fig.  8,  a  10  times  decrease  of  kc  and  a  5  times  in¬ 
crease  of  Zc  lower  the  cathode  potential  by  approximately  the 
same  amount.  Thus,  the  relative  effect  of  the  catalyst  layer 
thickness  is  higher  compared  to  the  protonic  conductivity. 

Lastly,  the  effect  of  membrane  ohmic  loss  as  described  by 
Eq.  ( 1 3)  is  analyzed  in  addition  to  the  cathode  polarization  be¬ 
haviour,  by  comparison  with  recently  obtained  experimental 
data.  The  polarization  curves  for  a  25  cm2  active  geometric 
area  fuel  cell  composed  of  Toray  carbon  paper  gas  diffusion 
layer,  0.5  mg  cm-2  Pt  load  (i.e.  20  wt.%  Pt  on  Vulcan  XC-72) 
and  Nafion  1 15  were  recorded  at  353  K  and  2.5  atm  pressure 
on  both  air  and  FU  side  under  two  conditions:  un-hydrated 
new  MEA  and  conditioned  (i.e.  fully  hydrated)  MEA.  The 
MEA  conditioning  was  achieved  by  mnning  the  cell  in  gal- 
vanostatic  mode  at  three  different  loads  until  the  cell  voltage 
drift  was  below  1  mV  per  hour  for  two  consecutive  hours  [48] . 

The  goal  was  to  evaluate  whether  the  developed  correlat¬ 
ing  equations  could  predict  the  single-cell  behaviour  under 
two  significantly  different  operating  conditions.  As  expected 
the  performance  of  the  unconditioned  MEA  was  poor,  for  ex¬ 
ample  at  7000  A  m-2  the  cell  voltage  was  about  0.3  V,  whilst 
the  fully  hydrated  MEA  at  the  same  current  density  gave 
0.57  V  (Fig.  9).  For  both  cases  the  modeling  equations  re¬ 
sulting  from  dimensional  analysis  showed  a  good  fit  with 
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Fig.  9.  Comparison  between  experimental  data  and  model  prediction  based 
on  correlating  equations  among  dimensionless  numbers  for  conditioned  and 
unconditioned  ME  A  operated  at  2.5  atm  air  and  H2  pressure,  353  K. 

the  experimental  data  (Fig.  9).  It  must  be  noted  that  in  the 
model  the  contribution  of  the  fL  oxidation  overpotential  was 
neglected  for  the  sake  of  simplicity. 

The  conditioned  MEA  results  (Fig.  9)  were  modeled 
by  coupling  the  membrane  ohmic  drop  given  by  Eq.  (16) 
with  V4  =  —  1  (i.e.  no  hydrostatic  pressure  difference  effect) 
with  the  following  catalyst  layer  polarization  models:  intrin¬ 
sic  kinetic  control  (Eq.  (30)  with  to  =  6  x  10-4  Am-2  and 
b  =  0.085  V)  up  to  750 Am-2,  and  mixed,  kinetic  and  mass 
transfer  control  expressed  by  Eq.  (32)  with  yio  =  2,  with¬ 
out  gas  diffuser  flooding  for  9500  >  i  >  750  A  m-2  and  with 
gas  diffuser  flooding  (Eq.  (44),  3/22  =  —0.65  and  yo,d  =  30)  for 
i  >  10,000  A  m-2.  The  abrupt  voltage  drop  in  the  case  of  the 
unconditioned  MEA  on  the  other  hand,  could  be  described 
by  the  previously  employed  membrane  ohmic  drop  equation 
( 1 6)  in  conjunction  with  the  regime  III  catalyst  layer  polariza¬ 
tion  equation  (38)  expressing  mixed  kinetic  and  ohmic  drop 
limitation  (kc  =4  S  m-1),  coupled  with  diffuser  flooding  for 
i  >  4000  A  m-2.  Thus,  it  could  be  concluded  that  the  correla¬ 
tions  among  dimensionless  numbers  could  lead  to  practically 
useful  models  of  the  MEA  performance  and  moreover,  they 
could  offer  insights  into  the  phenomena  responsible  for  the 
polarization  behaviour  without  increased  mathematical  com¬ 
plexity. 


the  catalyst  layer  under  combinations  of  electrode  kinet¬ 
ics,  oxygen  and/or  ionic  transport  limitations  (Eq.  (25)) 
and  the  mass  transfer  in  the  gas  diffuser  (Eqs.  (44)  and 
(45)).  The  dimensionless  correlations  were  validated  against 
experimental  data  and/or  independent  mathematical  model¬ 
ing  results  reported  in  the  literature. 

A  dimensional  analysis-based  approach  coupled  with  a 
mathematical  model  of  the  fuel  cell  systems  could  be  use¬ 
ful  to  quantify  and  analyze  new  experimental  data  for  a 
wide  range  of  conditions  by  determining  the  exponents 
and  coefficients  y,  in  the  corresponding  monomial  correla¬ 
tions.  It  is  proposed  that  using  representative  single  cells 
for  a  specific  fuel  cell  stack  development,  a  data  bank  of 
yi  exponents/coefficients  could  be  generated  for  a  variety 
of  operating  conditions,  by  fitting  the  experimental  data  to 
the  various  correlating  functions.  In  the  case  when  y,  is 
dependent  on  the  employed  range  of  the  variables  (e.g., 
low  current  density  experiments  yield  a  set  of  values  dif¬ 
ferent  than  those  obtained  at  high  current  densities),  in¬ 
terpolation  of  the  dimensionless  correlation  could  be  per¬ 
formed  according  to  the  method  developed  by  Churchill 
and  Usagi  [27].  The  monomial  correlations  for  the  MEA 
in  conjunction  with  the  experimentally  determined  y,-  expo¬ 
nent/coefficient  data  bank  could  be  easily  implemented  in 
more  complex,  design  specific,  mathematical  models  of  the 
fuel  cell  stack. 

Lastly,  while  the  present  dimensional  analysis  has  been  ap¬ 
plied  to  the  equations  of  the  Bernardi-Verbrugge  PEM  fuel 
cell  model,  it  is  envisaged  that  a  similar  procedure  could  be 
employed  for  other  models  of  the  PEM  fuel  cell  and  for  differ¬ 
ent  fuel  cell  types  such  as  solid  oxide  and  molten  carbonate. 
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5.  Conclusions 

The  Quraishi-Fahidy  dimensional  analysis  technique 
was  employed  to  derive  characteristic  dimensionless  num¬ 
bers  for  the  membrane-gas  diffusion  electrode  assembly 
in  polymer  electrolyte  fuel  cells.  It  was  shown  for  the 
first  time  that  monomial  correlating  functions  among  di¬ 
mensionless  ratios  could  be  used  to  describe  concisely 
and  in  an  easy  to  use  form,  the  non-linear  dependence 
of  the  membrane  ohmic  drop  as  a  function  of  operat¬ 
ing  conditions  (Eq.  (12)),  the  polarization  behaviour  of 


Appendix  A.  Membrane  characteristics 

Table  A1  lists  the  physical  properties  and  the  ionic  (fixed¬ 
negative  charge-site)  concentrations  for  the  1100  equivalent 
weight  series  Nafion  (i.e.,  N1 12,  NE1 135,  N1 15  and  N1 17) 
and  Dow  membranes  relevant  to  the  present  work. 

The  fixed  charge-site  concentration  Cf  and  the  equivalent 
weight  EW  are  given  by  the  following  equations: 


AC  BW 
Ufl  Zm 


(A.l) 
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and 


EW  =  — , 
AC 


(A. 2) 


where  AC  is  acid  (or  proton  exchange)  capacity  (eqkg-1), 
BW  the  membrane  basis  weight  (kg  m-2)  and  Zm  the  wet 
membrane  thickness  (m). 

The  hydraulic  and  electrokinetic  permeabilities  of  the 
membranes  at  353-368  K  were  [5]:  kp  =  1.8  x  10“ 18  m2  for 
bothNafion  and  Dow,  while  k&  =  7.19  x  10-20  m2  forNafion 
and  1.13  x  10“ 19  m2  for  Dow  membranes,  respectively. 

For  N117  the  membrane  specific  conductivities  under 
open-circuit  condition  kq  for  Figs.  2  and  3  were  calculated 
employing  X=  11  and  16,  respectively,  in  Eq.  (11).  Thus  ko 
is  equal  to  13.2  Sm-1  for  conditions  relevant  to  Fig.  3  and 
10.3  S  irf  1  for  Fig.  2,  respectively.  For  the  rest  of  the  Nafion 
membranes  kq  was  obtained  by  applying  the  thickness  cor¬ 
rection  factors  with  respect  to  ato.ni  17  derived  from  the  work 
of  Biichi  and  Scherer  [34]: 


K'o.Nin 

K0,N112 

K0.N117 

K(),N115 


1.26, 

1.07. 


K0.N117 

*3),NE1135 


=  1.18 


and 


(A. 3) 


In  the  case  of  the  Dow  membrane  at  368  K,  it  was  assumed 
that  kq  is  about  80%  higher  than  the  conductivity  of  N117 
under  the  same  conditions  [5]. 
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